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Abstract 

The library LINPRO which provides the solution to the linear inverse prob- 
lem for data contaminated by a statistical noise is presented. The library 
makes use of two methods: Maximum Entropy Method and Singular Value 
Decomposition. As an example it has been applied to perform an analytic 
continuation of the imaginary time propagator obtained within the Quantum 
Monte Carlo method. 
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1. Program Summary 

Title of the program: LINPRO vl.O 
Catalogue number: .... 

Program obtainable from: CPC Program Library, Queen's University of 
Belfast, N. Ireland (see application form in this issue) 
Licensing provisions: GNU Lesser General Public Licence. 
Distribution format: tar.gz 
Programming language: C++ 

Technical and API documentation: Yes, in HTML format 

Computer: LINPRO library should compile on any computing system that has 

C++ compiler. 



This paper and its associated computer program are available via the Computer 
Physics Communications homepage on ScienceDirect 
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Operating systems: LINUX or UNIX. 

Tested with compilers: GNU Compiler g++, Intel Compiler icpc. 

External libraries: OPT++: An Object-Oriented Nonlinear Optimization 

Library jsj (included into distribution). 

No. of lines in distributed program, source files only: 8 517. 

Nature of problem: LINPRO library solves linear inverse problem with an 

arbitrary kernel and arbitrary external constraints imposed on the solution. 

Solution method: LINPRO library implements two complementary methods: 

Maximum Entropy Method and SVD method. 

2. Linear inverse problem 

2.1. Formulation of the problem 

The inverse problem considered here is of the form: 

/oo 
Kix,y)A{x)dx, (1) 
-oo 

where y G [a, /3) and the kernel i^' is a known, real function, sufficiently regu- 
lar, although not necessarily smooth. The function G is known, and is repre- 
sented by a finite number Nr of values at a given set of points: {yi,y2, ...,yN^). 
The values G{yi) = Gi and G = (Gi, G2, ■ ■ ■ , Gjy^)'^ will be called the data 
and the data vector, respectively. These values are assumed to be in addi- 
tion affected by a noise of statistical origin and has to be treated merely as 
approximations of the true values. The unknown function A will be called 
the object, irrespective to its physical nature. The object is assumed to be 
nonzero only within a finite interval {a,b), although a and b are in general 
unknown. Moreover A may be a subject of additional constraints of the form: 

/oo 
gi{x)A{x)dx = Ci, i = 1,2,...,L (2) 
-00 

and 

A{xj) e[lj,Uj], j = l,2,...,M, (3) 

where functions g^ and values Cj are known, Ij and Uj indicate the lower and 
upper bound imposed on the object at some point xj. 
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2.2. Normal solution 

Since the function G is known for the finite set of argument values the 
hnear inverse problem ([1]) in practice reduces to its discretized counterpart: 

/OO /"OO 
K{x,y.i)A{x)dx= / K*{x)A{x)dx = iKi,A), (4) 
-OO J — OO 

where (■, ■) denotes the inner product. The object A can be treated as an 
element of A^- dimensional Hilbert space An (in general = oo). Note that 
due to discretization the kernel functions Ki{x) span only M-dimensional 
subspace Am {M < Nr) of the space An- It makes the inverse problem ill- 
posed, as there exists an infinite class of solutions satisfying Eq. (j4]). Indeed, 
let us expand Ki and A in an orthonormal basis {uk}^^i in An'- 

AI 

= fikUkjx), (5) 

k=l 

N M N 

k=l k=l k=M+l 

= Ap{x) + A^{x), (6) 

where Ap G Am represents the projection of A onto the M-dimensional 
subspace of An and A± is the remaining part, orthogonal to Ap: {Ap, A±) = 
0. Substituting the above expansions to Eq. one gets 

M N M 

'^i = ^Yl /ifc«/(wfe, ui) = f-^ak, (7) 

k=l 1=1 k=l 

where we have used the property {uk,ui) = 6ki- The last equality shows that 
Gi is independent of A±, since Gi = {Ki,A) = {Ki,Ap). It implies that the 
data vector G allows only for the reconstruction of Ap. The solution Ap 
with the minimal norm is a unique element of the subspace An and is called 
the normal solution 

In the case of data contaminated by a statistical noise the solution of the 
problem is also affected by uncertainties. Below we present two strategies 
which allow us to deal with such problems: 

1. the singular system analysis which uses the Singular Value Decompo- 
sition (SVD) to determine Ap and subsequently decrease uncertainties 
of the normal solution by incorporating constraints imposed on A, 
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2. the Maximum Entropy Method (MEM), which finds the most proba- 
ble solution, under the condition that data represent random numbers 
normally distributed around the true values, and that certain objects 
A are more probable than the others (so called a priori information 
about the object A). 

3. Singular system analysis 

3.1. SVD Method 

The normal solution Ap can be determined usingthe singular value de- 
composition of the integral kernel in Eq. (jl]) [H, l^l Let us rewrite it in 
the form: 

G = ICA. (8) 

The kernel functions Ki{x) span M-dimensional subspace Am and there- 
fore G has only M independent elements Thus G is an element of M- 
dimensional vector space Gm^ ■ /C is an integral operator which transforms 
an object from ^^r-space into a vector of the data space Gm^- The operator 
can be treated as a rectangular matrix of dimension A^^ x A^. We can define 
also a conjugate operator JC^ which transforms vectors from ^^^-space into 
^TV-space using the relation: 

{}Cu,v) = {u,IC^v), (9) 

where u G An, v G Gj^f and the inner product in the data space is defined 
as follows: 

{v,v') = j2^iv:, (10) 

1=1 

which is a useful definition in the case of uncorrelated dat80. The conjugate 
operator can be treated as a rectangular matrix of dimension N x Nt. 
Consequently the operator /C/C^ is represented by a square matrix of dimen- 
sion A^^ X Nr- Matrix elements of the operator /C/C^ are simply given by 

{}C}C% = {K,,K,). (11) 



^In the case of correlated data it is more appropriate to define the inner product as 
{v,v') — X]fj-=i ^iWijVj, where the matrix W is the inverse of the covariance matrix. 
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Performing the diagonalization of the matrix /C/C^ enables to determine the 
dimension of the subspace spanned by the kernel functions Ki. Indeed, the 
operator IClO has M positive eigenvalues {A^}^^, where M is the rank of the 
operator KK) . Corresponding eigenvectors {vi}fLi form a basis in the data 
space. Conjugate operator IOK, which acts in the object space An, among 
its eigenvalues has the same positive eigenvalues as the operator KK) and 
its eigenfunctions {wj}^! form the basis of ^M-space. The eigenvalues {Af }, 
the eigenvectors {vi\ and the eigenfunctions form a singular system of 
the operator /C satisfying the shifted eigenvalue problem: 

)Cui = XiVi, K)vi = XiUi. (12) 

The numbers Aj are singular values, and Wj, Vi are singular functions and 
singular vectors, respectively. The definition of the conjugate operator lO 
and equations (1121) allow to express the singular functions in the form: 



1 

= Kk{x){vi)k, (13) 



fc=i 



where {vi)k denotes k-th element of vector (vi). 

The singular system forms a suitable basis for expansion of the unknown 
object [l|: 

M 

Ap{x) = J2bMx), (14) 

i=l 

where the expansion coefficients are given by 

._iv.,G) ^^^^ 



3.2. Data with noise 

The solution given by Eq. f[T^ can be used only in the case of noiseless 
data. In the case when data vector G is known with some uncertainty AG 
the above algorithm becomes numerically ill-conditioned [H, Q, |3|. Note that 
the singular values {Aj}, the singular vectors {vi\ and therefore the singular 
functions {ui\ are known exactly since they are fully determined by the kernel 
functions Ki. Errors AG affect only the expansion coefficients, which will be 
the subject to some uncertainty A6j = (t?j, AG)/Aj. To perceive the origin 
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of numerical instabilities let us arrange the set of singular values {\i}fii in 
descending order: Ai ^ A2 ^ • • • ^ Am- Clearly, with a decreasing singular 
value the contribution of the statistical noise to is amplified: 

A, -^0 Ak= ^oo. (16) 

Ai 

Practically it means that in the object space Am there exist "directions" 
which are invisible for the SVD method, namely, the expansion coefficients 
cannot be determined with sufficient accuracy starting from some i index. 
Simultaneously the corresponding singular functions Ui become rapidly oscil- 
lating with an increasing index i (number of nodes of i-th singular function is 
i — 1) The functions associated with smaller values of Aj are responsible 
for reconstructing more subtle details of the solution. Since large uncertain- 
ties of coefficients in general yield to strong fluctuations of the solution, one 
of the standard methods is to remove all such strongly fluctuating terms and 
include only those for which bi are determined with satisfactory accuracy: 

Mcut 

Muti^) = ^biU,{x). (17) 
1=1 

This approach leads to the so called truncated SVD method (TSVD). In 
practice the truncation parameter is chosen in such a way to remove all terms 
for which the ratio (cut-off parameter) A^/Ai is smaller than 
ensures that the solution Ap^^^ reproduces data Gi within its error bars and 
prevents the inclusion of unjustifled structures into the solution j2|, [sj. 

3.3. Incorporating a priori information 

The reconstruction quality of the SVD method decreases significantly if 
data are affected by even a relatively weak noise. It turns out however that 
the incorporation of a priori information can improve the reconstruction 
process 0, Isf. There are two types of the prior information: information 
concerning the support of the solution (interval where the solution is nonzero) 
and external constraints. The first type of information leads to the following 
modification of the original problem: 

/+00 pb p+oo 

K*{x)A{x)dx= / K*{x)A{x)dx= / K*{x)S{x,a,b)A{x), 
-00 Ja J —00 

(18) 
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where object A is assumed to be non zero in the interval {a,b). S{x,a,b) 
denotes the support function defined as 

which modifies the kernel functions for the SVD method. This modification 
has two major consequences. First, the singular values Aj decrease faster as 
the size |6 — a| of the support function gets smaller. It subsequently leads to 
smaller values of Mcut and in general decreases the reconstruction ability of 
the method. Second, however, the singular functions Ui become limited to the 
interval x G (a, b) and their zeros are spaced more closely. This implies that a 
smaller number of singular functions are needed to get the same accuracy of 
reconstruction as before. It turns out that the latter consequence dominates 
and a properly chosen support function increases reconstruction quality jsj. 

Within the SVD method it is also possible to generate the solution which 
satisfies integral constraints (|2]). This can be done using the fact that each 
solution of the form 

M 

Ap{x,{bi})=J2b^U^{x), (20) 
i=l 

where 6j G (6j — A6j, bi + A6j) reproduces the data Gi within its error bars. 
Hence choosing an appropriate set of the expansion coefficients {bi}fi^ one 
can try to reproduce constraints (it is not always possible since the normal 
solution need not fulfill the same constraints as the true solution) In 
general, the expansion coefficients {bi}^!^ which agree with the constraints 
are not unique. To distinguish between various possibilities one can de- 
fine the cost functional C[Ap], which has to be minimized to find the best 
set of coefficients. As a cost functional one can use statistics with a 
similar form like in the maximum entropy method (see next section). An- 
other possibility is to choose the cost functional as the norm of the solution, 
C [Ap] = 1 1 v4p 1 1 = a/ {Ap, Ap). This choice is in agreement with the spirit of 
an SVD approach, where the normal solution is defined as the solution with 
the minimal norm. Summarizing, the problem of determining the unknown 
object satisfying external constraints has been reduced to the optimization 
problem: 

Ap{x)=mmC[Ap{x,{k})] (21) 
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with external constraints: 

Vz = l,2,...,M : bi-Ab,^k^bi + Abi, (22) 

/ + 00 
gj{x)Ap{x)dx = Cj. (23) 
-oo 

4. Maximum Entropy Method 

4.1. General considerations 

Let us distinguish between the exact values of the function G which are 
unknown and fulfill Eq. ([1]) and their known approximations contained in a 

vector: G = (Gi, G2, • • • , G]^^Y . These values can be treated as a particular 
realization of random variables, which are assumed to be uncorrelatecl^ and 
have a normal distribution around the exact values Gi with a variance af. 

The probability of obtaining the particular realization G under the condition 

that the exact values are given by G reads 

p(G|G)ocexp|^-i£(^^i-^j j, (24) 

and the values Gi depend on the function A according to the relation 
This equation is subsequently discretized in a chosen interval (a, b) and be- 
comes a linear transformation: 

N 

G, = Y,K,,A„ (25) 

where Kij = K{xj,yi)Ax is a rectangular matrix A^^ x A^, Ax = xj — xj-i 
and Aj = A{xj). Points Xj are uniformly distributed over the interval (a = 

Xi,b = XN)- 

The estimator for the quantity A = {Ai, A2, . . . , A^)"^ is defined as the 
one which maximizes the conditional probability p{A\G). This in turn can 
be expressed by flMj) using Bayes' theorem: 

p(A\d) = E!M>p(^, (26) 
P(G) 



^The extension of the method to the case of correlated data is straightforward, but 
require additional information in the form of covariance matrix. 
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where on the rhs the dependence on A is included in G through the relation 
(123]) . The probability p{A) is a priori probability and may contain additional 
information about A including constraints The maximization of this 
probability (so called likelihood function) leads in practice to the condition: 

-^PiAG) = 0, j=l,..,N, (27) 

In the case of N < Nr and p{A) = const the above condition reduces to the 
least square problem with the solution: A = {K^K)^^G^K, where ai = a = 
const is assumed. 

Here we are interested in the case when > Nr and an additional prior 
information is needed. It is specified through the entropy S{A), where p{A) oc 
exp{S{A)). The completely non-informative entropy is of the form: 

where a > is arbitrary. It favors the solution A = const. Usually we have 
additional information about the structure of A which allows us to specify a 
model of A. In such a case the relative entropy can be constructed: 

SiA\M) = -aJ2^^I^, (29) 
i=i l^j=i 

where M. = {Aii, Mn)^ , -Mi = -Mixi) is an assumed model for A and 
'^f^i Ai = X^ili -^i- The model M. has to fulfill the constraints imposed on 
A. In order to be able to construct the entropy in the above form, requires the 
assumption of nonnegativity of A. In order to avoid a complicated notation 
we assume also that both A and A4 are normalized: Xlili ^« ~ Xlili ~ ^■ 
Clearly the entropy is maximized in the case when A = A4, although note 
that S{A\M) ^ S{M\A). 

The prior information provides additional conditions for A and makes the 
maximization of the likelihood function a well defined process with a unique 
solution. Clearly now: 

p(l|G)ocexp(-l^[^i-^] -c^f^Alog^], (30) 
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and the maximum entropy method leads to the maximization of the above 
function with respect to A [6[. Note that still one has a freedom of choosing 
the constant a > 0. It governs the relative importance of the two terms in 
the above expression and larger a favors the model over the data. 

Another extension of the above formulation which will be considered in 
the next section admits the possibility of having a class of models }A[x\ /), 
where / is a set of parameters describing admissible degrees of freedom of 
the model and thus defining a set of admissible models. 

Method of solution 
The quantity which has to be minimized as a result of the MEM reads: 

2 




The task of minimizing the function of variables, where N in practice may 
be of the order of 10^~^ is rather hard. Therefore we apply here the procedure 
described in Ref. [3] which replaces the minimization of the many-variable 
function by an iterative process of consecutive least square problems. Let us 
assume that F{A^) = min and represents the solution of the problem. 
We expand F around A^ up to the second order: 



F{A^ + 6A) = F{A^,A) 



where 



= M-A°(^l-log-^ + ^(log^) 1, (33) 

A = ^ 

The above expansion implies the method of solving the problem. Namely, in 
the step n we minimize 

i?(^o(n)^ ^(n)) with respect to A^") at fixed v4°("). This 
is equivalent to the least square problem. Then we define a new 
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^A^"'~^^^ + {1—^)A^"'\ where ^ G (0, 1). Such aprocedure leads to a convergent 
solution providing ^ is sufficiently small. As a starting condition one takes 
A(o) = A(') =M. 

The additional constraints ([2]) can be included by considering the modified 
function: 

L N 

G{A',A) = F{A',A) + J2e,{c, -J29^,,A,)^ (34) 

i=i j=i 

where Qij = gi{xj)Ax, and 6i are positive parameters governing the "stiff- 
ness" of the constraints and thus responsible for the accuracy at which con- 
ditions ([2]) are fulfilled. 

In the MEM approach we have improved the method of finding the so- 
lution by constructing a sequence of minimizations with a gradually refined 
model. In this case the model is of the form Ai{x; f) and thus represents 
a class of models defined by parameters / = (/i, /g). At the end of each 
minimization process described above the result has been used to define a 
new model which maximize the overlap with respect to parameters /. 

The above quantity is clearly nonnegative and moreover 0{f) G [0, 1]. It is 
equal to unity if Ai = Aii- After the maximization of the overlap the new 
minimization process is started as described above. The procedure has been 
continued until the value of \A^"'^ — A^'^~^^\ < e, with an admissible tolerance 
e > 0. This strategy will be called as "self-consistent" Maximum Entropy 
Method. 



5. Structure of the library 

5.1. General overview 

LINPRO is an object-oriented library for solving linear inverse problems 
written in C++. As an optimization engine it uses OPT++ 2.4 library jsj. 
The aim of the LINPRO library is to collect in one place various algorithms for 
solving inverse problems and provide unified and user friendly programming 
interface to all of them. The library can be used to solve the problem with an 
arbitrary kernel defined by the user. For the MEM the package provides an 
interface for defining the arbitrary default model as well as a class of default 
models parametrized by a set of parameters. 
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5.2. Installation and technical documentation 

To install the library unpack tarball and follow instructions contained in 
INSTALL file. The distribution contains also folder doc where the technical 
and the API documentation is located. The documentation is generated in 
the HTML format, the master file is index.html. The user will find the codes 
which solve the artificial problem (as presented below), in the attached folder: 
examples. 

5.3. Inverse problem solvers 

Within LINPRO library algorithms for solving the inverse problem are 
called solvers, represented by InverseProblemSolver class. Fig. [1] presents 
the inheritance diagram of available solvers. 



GenericMEMSolver 



Ej^pMEMSolver 



MEMSolver 



InverseProblemSolver 



SCMEMSolver 



SVDSolver 



StdMEMSolver 



Figure 1: (Color online) Inheritance diagram for InverseProblemSolver. 



The solvers represents algorithms: 

• StdMEMSolver - standard implementation of the Maximum Entropy 
Method. To minimize the likelihood function the algorithm reduces 
the minimization problem to the iterative process of consecutive least 
square problems, as described in section 14.21 

• MEMSolver - implementation of the Maximum Entropy Method. It 
minimizes the likelihood function using the nonlinear interior-point 
method. 

• ExpMEMSolver - implementation of the Maximum Entropy Method, 
where object A{x) is parametrized by the formula A{x) = A4{x) exp f{x) 
where Ai{x) is a model function and f{x) is determined by the solver. 
Such a substitution is often used to eliminate the term \ogA{x)/A4{x), 
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which in specific situations is a source of optimizer instabihties (for ex- 
ample such instabihties can occur if a chosen model A4 (x) is very close 
to zero for some values of x, and due to finite precision is treated as 
zero). To minimize the likelihood function the nonlinear interior-point 
method is used. Since this solver is much slower than standard solvers 
it should be used in cases when StdMEMSolver and MEMSolver do not 
converge properly. 

• SVDSolver - implementation of the SVD method. 

• SCMEMSolver - implementation of the self-consistent engine for Maxi- 
mum Entropy Method. It works with each solver belonging to GenericMEMSolver 
branch. The scheme of the solver algorithm presents Fig. [21 



SCMEMSolver 



Use 

GenericMEMSolver 
to solve problem 
for fixed model 
A 



Update model to 
maximize overlap 
with solution 



Check 
self-consistency 
criteria 



STOP 



Figure 2: (Color online) The algorithm scheme for SCMEMSolver. 



6. Example of physical application 

6.1. Problem formulation 

As an example the package has been applied to extract the spectral weight 
function A{x) through the analytic continuation of the imaginary time prop- 
agator G{y): 

G(y) = -— r dxA{x) ^^P(-;^^) (36) 
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By definition, A{x) fulfills the following constraints: 



A{x) > 




(37) 




1 + exp{x(3) 



1 



G(/3). 



(38) 



This problem is frequently encountered in Quantum Monte Carlo simula- 
tions, which by construction produce data affected by the statistical noise 
[il, [lO]- In order to check the reconstruction ability of the package an artifi- 
cially generated data for the imaginary time propagator has been used. The 
application of the package to the real physical data can be found in Refs. 



The artificial spectral function has been chosen in the form: 



where N{x; fi, a) is the normal distribution function with the mean /i and the 
standard deviation a. Subsequently the values of the imaginary time prop- 
agator has been generated using the relation (l36l) for uniformly spaced 
data points in the interval [0,/3 = 10]. 

6.2. Data without noise 

As a first step the quality of reconstruction of the object as a function of 
the number of data points N^- has been tested. It was found that in order 
to reproduce the original object with a satisfactory accuracy, one has to use 
Nr ^ 20 data points; see Fig. [3l Further increase in the number of data 
points does not improve significantly the conformity between the solution 
and the object. It is related to the fact that the dimension of the subspace 
Am increases linearly with an increase in the number of data points up to 
Nr = 20 and then it saturates (for the presented example the dimension 
of ^Af-space is 20 for A^^ = 20 and 25 for A^^ = 100). Note also that 
the projected solution produced by the SVD method provides a very good 
approximation of the "true" object. 

6.3. Data with noise 

In order to test the ability of reconstruction in the presence of noise each 
value Gi has been perturbed by the Gaussian noise of zero mean value and 



A{x) = ^N{x; -1.5, 0.5) + ^N{x; 2.0, 0.7) 



(39) 
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Figure 3: (Color online) The reconstruction of the artificial object function A{x) by the 
SVD and MEM methods. The reconstruction is performed using Nr uniformly spaced 
data points d in the interval [0, /3 — 10] (noiseless data). A very good agreement between 
the normal solution and the original object is achieved if Nr ^ 20. The root mean square 
error (RMSE) for both methods is displayed at the bottom of the figure. 

the standard deviation equal to 1% of Gi. In the case of an SVD method the 
object has been reconstructed using Eq. ( IT71) for various cut-off parameters 
Ai/Ai (Ai is the highest singular value); see Fig. HI Note that solutions 
with a cut-off parameter bigger than the relative error 0.01 do not guarantee 
the reproduction of the imaginary time correlator within its error bars. It 
is clearly seen that with a decreasing value of the cut-off parameter the 
quality of reconstruction increases, and the "optimal" cut-off parameter is 
Aj/Ai ~ 0.01. Further decrease in the cut-off leads to the inclusion of an 
unjustified structure (strong fluctuations) into the shape of the reconstructed 
object. 

In the case of MEM the quality of reconstruction is a function of a pa- 
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Figure 4: (Color online) Reconstruction of the artificial object function A{x) by TSVD 
method. Reconstruction is performed using A^,- = 25 uniformly spaced data points Gj 
in the range [0,/3 = 10], perturbed by Gaussian noise AGi — A/'(0, Gi/100). In the 
reconstruction procedure only those terms were included for which Xt/Xi is larger than 
a given cut-off parameter. Mcut denotes the number of the singular function included in 
TSVD expansion. 



rameter; see Fig. [51 A class of assumed models has been chosen according to 
the prescription: 

A^(x;ci,C2,/ii,/i2,cri,(T2) = CiN{x; fii,cri) + C2N{x; ^12,(^3)- (40) 

It was found that there exists a critical value of a parameter which sep- 
arates "smooth" and "rigged" solutions. It corresponds to the value which 
minimize the total MEM errors as discussed in Ref. ll|. Moreover, it turns 
out that the "self-consistent" algorithm always converge to the same solution 
irrespective of initial values of parameters which define the class of models. 
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Figure 5: (Color online) Reconstruction of the artificial object function A{x) by self- 
consistent MEM for different values of parameter a. 

6.4- Impact of external constraints and a priori information 

In the following the impact of the external constraints on the reconstruc- 
tion quality has been tested for both MEM and SVD methods; see Fig. El 
The external constraints influence strongly the SVD method. Note also that 
the solution produced by an SVD method is not an accurate reconstruction 
of the input spectral function. It is due to the fact that the solution produced 
by the SVD method is a projection of the "true" spectral function onto the 
"visible" subspace, where the problem is well posed. The main advantage of 
an SVD approach is that it does not require any a priori information. Conse- 
quently the SVD solution can deliver very useful information concerning the 
default model or a class of default models for MEMs. In this particular test 
the SVD solution suggests that it is profitable to choose the default model 
as a combination of two Gaussians (left panel), given by Eq. (HOj) . One can 
also use the SVD solution as a default model for MEMs. 
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The tests presented above suggest that the maximum entropy method 
combined with "self-consistent" engine provides the most accurate solutions. 
Even if the class of the models is not correctly prepared, the self-consistent 
solution still well reproduces the input object. The right panel presents the 
case where the class of default model Gaussian functions N{x; n, a) was used. 
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Figure 6: (Color online) The reconstruction ability of the spectral function for the full 
problem (data with noise -|- external constraints) of the SVD and MEM methods. The left 
panel shows the solution of the self-consistent MEM with a combination of two Gaussians 
functions as a default model class. The right panel shows the solution of the self-consistent 
MEM with Gaussian functions as a default model class. 

Therefore the best methodology of producing the solution is suggested to 
be the following: 

1. Create an SVD solution and apply it to construct the class of default 
models Ai{x] /); 

2. Use the "self-consistent" MEM with constructed class of models A^(a;; /) 
to produce final solution. 

6.5. Resolution limit 

In the case of physical applications where the object is associated with 
the spectral weight function, an extremely important question needs to be 
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Figure 7: (Color online) Left panel: the resolution limit as a ftinction of temperature 
T = 1/13 (red solid line). Above the limit it is possible to reconstruct the gap (defined as a 
distance between two peaks). The sketch of a typical evolution of the physical gap (dotted 
blue line) is also plotted. The gap can be properly reconstructed up to T* temperature. 
T* denotes the temperature for which the gap vanishes, assuming that the reconstruction 
provide an exact object. Right panel: the value of the gap Aj- reconstructed by the 
SVD method versus the true value Aq. The support function is ^(a;, — 2,2). Arrows 
indicate the minimal value of the gap (resolution limit), where the bimodal structure of 
the reconstructed spectral function appears. 

answered: is the spectral function unimodal or bimodal? It is well known 
that the distinct peaks of a bimodal spectral function may be overlooked 
during the reconstruction process if the distance between peaks is smaller 
than some critical value, which defines the resolution limit. 

To quantitatively estimate the resolution limit, the artificial object func- 
tion consisting of two delta functions separated by 2Aq distance have been 
considered. Namely, A (x) = S{x+Aq)+6{x—Aq). For this function the imag- 
inary time propagator has been generated for A^^- = 25 uniformly distributed 
data points in the interval [0,/3], where now /3 is treated as a parameter (in 
this case the inverse of /3 has the physical meaning of temperature). The 
resolution limit Aq™'"^ is defined as a minimal value of Aq for which the bi- 
modal structure of the object can still be reconstructed and in general is a 
function of (3. 

Results of the presented tests are shown in Fig. [71 For both methods 
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(SVD and MEM) the existence of the finite resolution Umit has been found. 
It increases with a decreasing temperature T = (left panel). Within the 
presented approach it is possible to reconstruct the gap (defined as a distance 
between two peaks) only for the temperatures for which it is larger than the 
reconstruction limit. 

Let us consider a process of reconstructing the physical gap, which is a 
decreasing function of temperature and eventually vanishes at some temper- 
ature T*. At a certain temperature T* < T* the gap becomes comparable 
with the reconstruction limit. Up to this temperature the reconstructed value 
of the gap Ar agrees very well with the true value Aq (right panel). At the 
temperature T* the value Ap drops to zero. It means that the methods pro- 
vide in practice a lower bound for the temperature at which the true gap 
vanishes. 

7. Conclusions 

Library LINPRO for solving arbitrary linear inverse problems with exter- 
nal constraints has been presented. The library uses the Maximum Entropy 
Method and the SVD methods. An object-oriented implementation ensures 
that the package acquires a unified and user friendly interface. As an example 
we have applied the new package to solve the typical problem of computa- 
tional physics: analytic continuation of imaginary time propagator to real 
frequencies. 
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